% phi_res in data_from_stata2
% res1 and res2  in bias_comp

clear all

nreps = 1000;

bootgamma = zeros(5 ,nreps);
bootphi =  zeros(10 ,nreps);
bootres1 =  zeros(2 ,nreps);
bootres2 = zeros(10 ,nreps);
bootimpb = zeros(5 ,nreps);
bootrobres = zeros(9 ,5,nreps);

save bres bootphi bootres1 bootres2 bootimpb bootgamma bootrobres

for b = 1:nreps
b
%set seed
seed = 45454 + 10*b;
rand('state', seed);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  load hedonic data
load hed_data_log

dumyear = dummyvar(year);
 
M = max(hed_m);
b1_mm_l = zeros(M,1);  b2_mm_l = zeros(M,1);
vMin = zeros(M,3); vMax = zeros(M,3);
meanxb = zeros(M,1);

load datafromstata residp Nm
for m = 1:M
    index = find(hed_m==m);
    logp(index) = logp(index) +  randsample(residp(index),Nm(m),'true');
end


hedrestable = zeros(7,M);
hedsetable = zeros(7,M);
Nm = zeros(M,1);
Rsq = zeros(M,1);


for m = 1:M
    index = find(hed_m==m);
    temptract = tract(index);
    dumtract = dummyvar(temptract);
    temp = mean(dumtract);
    index2  = find(temp>0);
    dumtract = dumtract(:,index2);    
    dumtract = dumtract(:,2:end);
    cons = ones(length(index),1);
    [temp ,~,r,~,stats] = regress(logp(index), [  v_crime(index) hv_crime_sq(index) X(index,:) cons dumtract dumyear(index,2:end) ]); %dumyear(index,:)
    se = sqrt(diag(inv(length(cons))*(r'*r)*inv([  v_crime(index) hv_crime_sq(index) X(index,:) cons dumtract dumyear(index,2:end) ]'*[  v_crime(index) hv_crime_sq(index) X(index,:) cons dumtract dumyear(index,2:end)])));
    hedrestable(:,m) = temp(1:7);
    hedsetable(:,m) = se(1:7);
    Nm(m) = length(r);
    Rsq(m) = stats(1);
    b1_mm_l(m) =temp(1);
    b2_mm_l(m) =temp(2);
    meanxb(m) = mean([X(index,:) cons dumtract  dumyear(index,2:end)])*temp(3:end);
    
    vMin(m,1) = min(v_crime(index));  vMin(m,2) = prctile(v_crime(index),1); vMin(m,3) = prctile(v_crime(index),5);
    vMax(m,1) = max(v_crime(index));  vMax(m,2) = prctile(v_crime(index),99); vMax(m,3) = prctile(v_crime(index),95);
end
fvMin = vMin;
FvMax = vMax;
% Xvars are age lotsize sqft num_rooms
matrix = zeros(12,5);
matrix(1,:) =  1000*hedrestable(1,:);  matrix(2,:) =  1000*hedsetable(1,:); %vc
matrix(3,:) =  1000000*.5*hedrestable(2,:);  matrix(4,:) =  1000000*.5*hedsetable(2,:); % vcsq
matrix(5,:) =  1000*hedrestable(5,:);  matrix(6,:) =  1000*hedsetable(5,:); % sqft
matrix(7,:) =  1000*hedrestable(4,:);  matrix(8,:) =  1000*hedsetable(4,:); % lot
matrix(9,:) =  hedrestable(3,:);  matrix(10,:) = hedsetable(3,:); % age
matrix(11,:) =  hedrestable(6,:); matrix(12,:) = hedsetable(6,:); % num rooms




% matrix = matrix;
% 
% rowLabels = {'Violent Crime Rate', '', 'Violent Crime Rate Squared', '', 'House Square Footage', '', 'Lot Square Footage', '', 'House Age', '', 'Number of Rooms', ''};
% matrix2latex_a2(matrix, 'hedres.tex', 'rowLabels', rowLabels, 'columnLabels', [], 'alignment', 'c', 'format', '%-6.5f', 'size', []);


clear  index temptract dumtract temp index2 bint r rint stats v_crime hv_crime_sq logp year tract X




%  load transition prob data
load tprob_data

M = max(tprob_m);
pxest = zeros(M,3);sigrhoest=zeros(M,1);
pxse = zeros(M,3);

load datafromstata residt Nmtprob
for m = 1:M
    index = find(tprob_m==m);
    v_crime(index) = v_crime(index) +  randsample(residt(index),Nmtprob(m),'true');
end


hedrestable = zeros(3,M);
hedsetable = zeros(3,M);
Nmtprob = zeros(M,1);
Rsqtprob = zeros(M,1);

vc = zeros(M,3,3);


for m = 1:M
    m;
    index = find(tprob_m==m);
    cons = ones(length(index),1);
    [temp ,~,r,~,stats] = regress(v_crime(index), [cons Lv_crime(index) tprob_t(index)]);
    pxest(m,:) = temp;
    sigrhoest(m) = sqrt(inv(length(index))*(r'*r));
    pxse(m,:) = sqrt(diag(inv(length(cons))*(r'*r)*inv([ones(length(index),1) Lv_crime(index) tprob_t(index)]'*[ones(length(index),1) Lv_crime(index) tprob_t(index)])));
    vc(m,:,:) = inv(length(cons))*(r'*r)*inv([ones(length(index),1) Lv_crime(index) tprob_t(index)]'*[ones(length(index),1) Lv_crime(index) tprob_t(index)]);
    hedrestable(:,m) = temp;
    hedsetable(:,m) = pxse(m,:);
    Nmtprob(m) = length(r);
    Rsqtprob(m) = stats(1);
end
pxest = pxest';
pxse = pxse';



% Xvars are cons lagged crime t
% matrix = zeros(6,M);
% matrix(1,:) =  hedrestable(1,:);  matrix(2,:) =  hedsetable(1,:); %cons
% matrix(3,:) =  hedrestable(2,:);  matrix(4,:) =  hedsetable(2,:); % lagged crime
% matrix(5,:) =  hedrestable(3,:);  matrix(6,:) =  hedsetable(3,:); % t
% 
% matrix = matrix;
% 
% rowLabels = {'Constant', '', 'Lagged Violent Crime Rate', '', 't', ''};
% matrix2latex_a2(matrix, 'tprobres.tex', 'rowLabels', rowLabels, 'columnLabels', [], 'alignment', 'c', 'format', '%-6.4f', 'size', []);




yvals = [1:1:19]'; x= 0; 
% with AR(1), xbar is linear in x. x can be evaluated at any value.
% if xbar is non-linear in x. x needs to be evaluated at every value.
bigtheta0 = zeros(5,length(yvals));
for y=1:length(yvals);
gamma = ones(M,1);
dgammadrho = zeros(M,1); %for delta method
beta = 0.95;
v_crime_bar = x*ones(M,1);
extp=x*ones(M,1);
for m=1:M
B = 1;
for t=2:7 
extp(m) =  pxest(1,m) + pxest(2,m)*extp(m)  + pxest(3,m)*(yvals(y)+t-1);
temp = beta^(t-1);
temp2 = (beta^(t-1))*( pxest(2,m)^(t-1));
temp3 = (beta^(t-1))*( (pxest(2,m)^(t-2))*(t-1));
B = B+temp;
v_crime_bar(m) = v_crime_bar(m) + temp*extp(m);
gamma(m) = gamma(m) + temp2;
dgammadrho(m) = dgammadrho(m) + temp3;
end
end
v_crime_bar = v_crime_bar/B;
gamma = gamma/B;
dgammadrho = dgammadrho/B;
theta0=v_crime_bar-x*gamma;
bigtheta0(:,y)= theta0;
end

theta0 = bigtheta0(:,10); 

aaa=diff(bigtheta0')';
theta0con = aaa(:,1);

%%%
phi_res(1,:) = bigtheta0(:,1)';  
phi_res(2,:) = theta0con';



gammase = hedsetable(2,:)'.*dgammadrho;  %delta method.  se(gamma(rho))  = se(rho)*(d gamma/ d rho)
clear  temp  r   v_crime_bar y x temp2 index extp m



%  load consump data
load consump_data

M = max(con_m);

for m = 1:M
    index = find( con_m==m & (b1_mm_l(m) +  b2_mm_l(m)*a < 0) );
    vMin(m,1) = min(a(index));  vMin(m,2) = prctile(a(index),1); vMin(m,3) = prctile(a(index),5);
    vMax(m,1) = max(a(index));  vMax(m,2) = prctile(a(index),99); vMax(m,3) = prctile(a(index),95);
end

Na = length(a);
clear consump_data

save datafromstata

rentd=zeros(Na,M); 
for m = 1:M
    rentd(:,m) = (b1_mm_l(m) + b2_mm_l(m)*a).*price;
end
temp = dummyvar(con_m);
rentd = sum(rentd.*temp,2);
rentd = .075*rentd;
index = find(rentd<0);
con_m = con_m(index);
rentd=rentd(index);
t = t(index);
save est_res_s2 rentd  M con_m t



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
bias_comp
Robustness_checks_b



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load datafromstata gamma phi_res 
load bias_comp res1 res2 imp_bias
load rob_res rob_res

load bres

bootgamma(:,b) = gamma';
bootphi(:,b) =  [phi_res(1,:)'; phi_res(2,:)'];

bootres1(:,b) =  res1';
bootres2(:,b) = [res2(:,1); res2(:,2)];
bootimpb(:,b) = imp_bias;

bootrobres(:,:,b) = rob_res;


save bres bootgamma bootphi bootres1 bootres2 bootimpb  bootrobres 

end

sephi = std(bootphi'); sephi = [sephi(1:5);sephi(6:10)]
seres1 =  std(bootres1')'
seres2 = std(bootres2'); seres2 = [seres2(1:5);seres2(6:10)]
seimpb =  std(bootimpb')
segamma =  std(bootgamma')

serobres = std(bootrobres,[],3)

save bres bootgamma bootphi bootres1 bootres2 bootimpb segamma sephi seres1 seres2 seimpb  serobres

%rerun estimation to finish with saved estimation results
clear
estimate
bias_comp
Robustness_checks



clear
load bres segamma sephi seres1 seres2 seimpb  serobres
load datafromstata gamma phi_res 
load bias_comp res1 res2 imp_bias
load rob_res rob_res
save vc_results
